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Abstract. We study the nonequilibrium dynamics of the linear to zigzag 
structural phase transition exhibited by an ion chain confined in a trap with 
periodic boundary conditions. The transition is driven by reducing the transverse 
confinement at a finite quench rate, which can be accurately controlled. This 
results in the formation of zigzag domains oriented along different transverse 
planes. The twists between different domains can be stabilized by the topology of 
the trap and under laser cooling the system relaxes to a helical chain with possibly 
nonzero winding number. Molecular dynamics simulations are used to obtain a 
large sample of possible trajectories for different quench rates. The scaling of the 
average winding number with different quench rates is compared to the prediction 
of the Kibble-Zurek theory, and a good quantitative agreement is found. 
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1. Introduction 

Gases of singly charged ions can be spatially confined by Paul or Penning traps 
[1]. When they are laser cooled to sufficiently low temperatures, they form Wigner 
crystals [5], which are examples of a self-organized systems. Wigner crystals have 
been observed and characterized in several atomic physics experiments [3] and have 
recently attracted considerable attention as a platform for investigating nonlinear 
dynamics close to criticality. Some examples of the studies of the nonlinear dynamics 
of ion crystals include the simulation of linear and nonlinear Klein-Gordon fields on a 
lattice [1] , the study of nucleation of topological defects [SJ |S] , as well as proposals to 
realize models related to dry friction and energy transport [7] and synchronization f5] . 
The idea behind these works is to use ion traps as a flexible platform for designing 
experiments aimed at realizing various models of nonlinear dynamics. The high degree 
of isolation of the ion crystals from the surrounding environment, which is usually 
achieved in experiments, implies also the possibility to enter the regime where quantum 
mechanical effects must be accounted for to describe critical phenomena [H IHl UHl E] > 
and where the quantum motion can be utilized for quantum information using trapped 
ions [13 [13]. 

In this paper we consider the non-equilibrium statistical mechanics of a chain of 
ions, when the frequency of the transverse potential is quenched across the linear to 
zigzag transition. This analysis has been previously performed by some of us in two 
dimension [5J |S] . In this paper we consider a trap which is invariant by rotation about 
the chain axis. Before we continue, let us provide some relevant information about 
the system we consider. In [M] it has been shown that the linear-zigzag transition is 
a continuous phase transition and it can be described bv means of Landau model. In 



Figure 1: Three snapshots of the central 200 ions in a 400 ion chain in a trap with 
periodic boundary conditions undergoing a phase transition from a linear to a zigzag 
configuration. The fraction of the chain displayed is mapped to an open chain for 
clarity. Under high transverse confinement, above the critical point, the chain is in a 
linear configuration (left). Whenever the transverse trapping frequency is decreased 
below the critical value, the chain develops a number of twists (middle). At the later 
stages of the evolution the twist and anti-twist annihilate and the ion chain relaxes to 
a configuration in which the total winding number remains constant (right). 
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a potential like the one here considered the symmetry breaking is associated with the 
formation of two chains on a plane, forming a polar angle which is distributed with 
equal probability in the interval [0,27r]. These considerations are strictly valid in the 
thermodynamic limit. When the number of ions is finite, the boundary conditions 
play an important role. In particular, for periodic boundary conditions the ground 
state in the zigzag phase may exhibit a twist. The assumption of periodic boundary 
conditions is consistent with experimental realizations, such as in storage rings [15(ll6j 
or in ion rings which can be observed in multipole traps [17l [18] . In these system, 
thus, the dynamics we discuss in this paper could be observed. 

The topology of the ion ring trap plays a crucial role in our study. If the transition 
is driven instead in a linear Paul trap under transverse isotropic confinement, zig-zag 
domains are formed in an early stage. Subsequently a local relaxation occurs in which 
the different domains coaligne and merge, resulting in the growth of a single domain 
that eventually extends over the whole chain. As a result, no information concerning 
the dynamics of the transition remains. 

Given an initial state of the system, which is a linear chain with some small 
random perturbation of the ion positions and velocities, it is usually considered 
impossible to predict the winding number of the final helical configuration without 
numerically solving the equations of motion, even if the dynamics is fully deterministic. 
However, it may be easier to make predictions about the expected average of the 
winding numbers resulting from many experimental runs. Hence, we shall pay 
particular attention to the universal features of the dependence of the average winding 
number of the helix on the quench rate of the phase transition. In annular Josephson 
tunnel junctions undergoing a thermal quench, the formation of fluxons (supercurrents 
carrying a single quantum of magnetic flux associated with a non-zero winding number 
of the order parameter) has attracted much attention, and was indeed experimentally 
addressed in [19]. In a recent work [20] the Kibble-Zurek (KZ) theory has been 
similarly applied to quantify the spontaneous circulation in a newborn Bose-Einstein 
condensate (BEG) created in a toroidal trap by quenching the chemical potential of 
the system at a finite rate. We shall see that ion chains constitute an experimentally 
advantageous platform to study the dynamics of spontaneous symmetry breaking, 
where KZ theory can analogously be used to predict the scaling of the winding 
number in the helical ion chain. We have used numerical simulations to verify the 
KZ prediction and we find that within the statistical uncertainty the KZ prediction 
holds. 

The paper is organized as follows. In section 2 we review the Ginzburg- Landau 
description of the linear to zigzag structural phase transition in three-dimensions. 
In section 3 we show that Kibble-Zurek mechanism (KZM) can be used to relate the 
average winding number to the quench rate. The results, which have been numerically 
evaluated, are presented in section 4 and compared with the predictions of the KZM. 
In section 5 the conclusions are drawn. 

2. Ginzburg-Landau model 

In |14j it was shown that near the critical point of the linear to zigzag phase transition 
the classical ground state and the low energy excitations of a chain of ions can be 
described using Ginzburg-Landau theory. Its time-dependent version governs the 
nucleation of kinks in the zig-zag phase [5j [6] . In this section, we present the time- 
dependent Ginzburg-Landau model, emphasizing the 3D nature of the transition, and 
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discuss some of its properties. The derivation starts by writing the potential energy of 
the system. Consider a system of ions confined to a space with periodic boundary in 
the X direction and experiencing harmonic potential in the transverse direction with 
frequency i/t- For N ions having mass m, charge e and coordinates = {xj,yj,Zj) 
the potential energy reads 



1 



.7 = 1 



2 47reo 



EE 

k j=ik 



(1) 



The ground state of the system is a linear chain, provided that is above a 
critical frequency, i^^'^'* = •\/7C(3)/2aJo — 2.051ujq; C(z) is the Riemann zeta function 
and the characteristic frequency is ujq = ^ j Aire^raa^ , where a is the spacing between 

the ions [M]. At frequencies just below the linear chain configuration is unstable, 
with the zigzag becoming a favourable configuration. In the rest of this section, a 
continuum description of the system will allow us to establish a close correspondence 
with the Ginzburg-Landau model. Nonetheless, we shall come back to the full model 
in the rest of the paper. 

In the vicinity of the critical point the energy of the system can be written as [14] 
y(2) + y(4) 



V 



2 a 



(2) 



g m 
1~a 



dx \5{y' + z') + h' {d^yY + h' {d.,zY 
dx {y^ + 2z2y2 ^ 



where 6 



,(c)2 



, and the formulas for h and g are given in |14) . but their exact 



values are not needed in the current work. The continuum approximation is valid if 
the spatial variation of the transverse displacement of the ions is small in relation to 
the ion spacing. 

We can define A = y + iz and rewrite in a more compact form 



V=-- 
2 a 



dx 



S\A\^ + h^d^Ad^A* 



(3) 



Figure [5] illustrates graphically the meaning of the order parameter A{x). \A{a) \ is the 
distance between the jth ion with the coordinate x — ja and the cc-axis. 9{x) = arg(A) 
is the angle between the xy plane and the plane containing the x-axis and the point 
where the jth ion is located. 

The equation of motion for the order parameter, A(x), can be obtained from (|3]) 
using the Euler-Lagrange equation and is given by 



d^A 



dt"^ ' ^ dt dx'^ 

where we have also added a friction force, with a friction coefficient 77, and a stochastic 
force e{x,t). The friction force and the stochastic force come from the interaction of 
the ions with a cooling laser (for further details check reference p| l2T|). 

Equation Q is known in the literature as the Ginzburg-Landau equation. It is 
extensively studied as it is one of the most important equations used to describe 
complex pattern formation. Along with the more general complex GL equation. 



dA 
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Figure 2: Diagram of the helical crystal along the a;- axis with only one ion at point 
P made visible. 9 = aig{A) and OP = \A\. Below the critical point the ion effectively 
moves in a wine bottle or a Mexican hat potential, which is also illustrated on the 
diagram. Above the critical point, the ion chain adopts a linear configuration in which 
the average value \A\ vanishes and 9 is undefined. 

this equation is used to describe, for example, phenomena such as nonlinear waves, 
second-order phase transitions, superconductivity and Bose-Einstein condensation 
P^ . Equation describes the second order phase transition at 5 = 0. Substituting 
the plane wave solutions, A{x,t) — ae^*^*^*'^^, into the linearized equation ^ gives 
the dispersion relation 



For (5 > 0, ImJ7 < and Reft = for all k i.e. all modes are decaying and the state of 
the system is stable. For 5 < there are k values for which Imfi > 0; these modes are 
unstable and will grow. The point 5 = is a critical point of as supercritical pitchfork 
bifurcation. 

Equation (j3|) arises naturally in physical systems where the order parameter has 
rotational symmetry. For example near the critical point, symmetry breaking in 
Josephson tunnel junctions [19] as well as Bose-Einstein condensation [23] (within a 
Gross-Pitaevskii description) have the same dimensionality and symmetry properties 
as the linear to zigzag transition in ion trap. As a matter of fact, all these systems 
belong to the mean- field theory universality class. 

Equation (|4]) allows for phase winding solutions, and crossing the structural 
transition leads to U{1) symmetry breaking, after which the ground state manifold 
is not simply connected, i.e. the first homotopy group of the manifold of degenerate 
ground states is 7ri(t/(l)/{l}) = 7ri(S'^) = Z, as in the formation of cosmic strings 
in the Abelian Higgs model [24] and the examples above. With periodic boundary 
conditions the total phase must be equal to 2ttW, where W is an integer known as 
the winding number, i.e. 



The stability of the phase winding solutions was analyzed in a number of works, 
see for example [IS] [25] . It is found that for each value of S there are phase winding 
solutions that are unstable because their winding number is either too high or too low. 



17 = ^ ± ^^/-7j + 4{h^k^ + 5). 
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When the number is too small, additional twist will form and when its too big then the 
twists will merge into one. The phase must be continuous and well defined whenever 
1^1 is not zero, so the winding number can only change discontinuously at some point 
in time and space where \A\ = 0. The scenario resembles the formation of fluxons 
in annular Josephson tunnel junctions under a thermal quench, where the winding 
number in the experiments reported so far is severely restricted to W = ±1 (trapping 
a single fluxon or anti-fluxon) [19] . It is also analogous to the creation of spontaneous 
quantized circulation in an atomic cloud undergoing Bose-Einstein condensation in a 
ring trap under a thermal quench [571 HO] • 

3. Formation dynamics and the Kibble-Zurek mechanism 

The aim of KZ theory [351 122] is to estimate the size of the domains that are formed 
when the the control parameter is quenched across the critical value. Knowing the 
average size of the domains after the phase transition allows to estimate the number 
of defects (e.g. domain walls). For the case of formation of helical structures it is also 
possible to use KZ theory to estimate the winding number |29j and this is the subject 
of this section. 

Suppose that the transverse confinement is changed in such a way that the quench 
of the control parameter is linear in time with amplitude Sq around the critical point 
t G [—tq/2,tq/2], which is reached a.t t — 



where the ground state of the ion chain is linear for S(t < 0), and a zig-zag configuration 
with = for S{t > 0) < 0. If the quench is very slow i.e. tq — >• oo, then 
the dynamics is adiabatic and we would expect to arrive at a state with the lowest 
possible W. This is because the first mode that becomes unstable during the quench 
through the critical point is the mode with the smallest wave number, fc, (or highest 
wavelength) and, since the quench rate is very slow with respect to the spectral gap, 
other modes are not be excited. For a fast quench, tq — 0, the higher modes become 
unstable at a rate faster than the time it takes for the system to relax to the lowest 
mode configuration. This results in helices with a higher winding number. It is known 
that in the vicinity of the critical point (S ~ 0) the correlation length and relaxation 
time have power law divergences 



in terms of the dimcnsionless control parameter e — 5 /{vl ^ and where the exponents 

V and z are the critical exponents determined by the universality class to which the 
system belongs. The critical exponents for the linear to zigzag phase transition are 

V — 1/2 and z — 2, which are those of a mean-field theory [14]. The KZ mechanism 
splits the dynamics of the system into three stages: adiabatic when e ^ 1, frozen 
around the critical point e = 0, and adiabatic again once e 3> 1. Within this impulse 
approximation the average size of the domains in the broken symmetry phase is given 
by the correlation length at the boundary between the adiabatic and impulse stages. 
The time scale, i, at which adiabaticity breaks down is the so-called freeze-out time 
and it can be found by matching the relaxation time to the time remaining until the 



m 
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crossing of the critical point as imposed by the quench i ~ ||(_( ■ This reads 

i^[r„r^)^. (10) 
This sets the average size of domains to 

^-^o[^)^. (11) 

In an ion chain with periodic boundary conditions the inter-ion spacing is homogeneous 
which leads to a constant critical frequency of the transverse confinement. As a result, 
the transition is crossed simultaneously all over the chain. Otherwise, the argument 
above is to be modified to take into account the inhomogeneous nature of the ion chain 
[SIB], see as weh [201 ISI]- 

The domains can be thought of as regions where the angle, 9, is approximately 
constant. Each domain will pick at random, therefore there will be about C/^ 
random variables, where C is the length of the ion chain. Let us assume that the 
random variables are uniformly distributed between —tt and vr and thus have a mean 
of zero and a variance of vr^/S. The distribution of the winding number, W, is then a 
convolution of C/^ uniform random variables. For large C, the central limit theorem is 
valid and W will have a Gaussian distribution with mean zero and variance Ctt^/ (3^). 
Hence, the root mean square of the winding number inherits a dependence on the 
quench rate from equation (fTT|) . 

and a{W) = y^ir /2{\W\) . For the ion chain, the divergence of the correlation length 
derived from the two-point correlation function takes the form that ^ = aujQ / ^/S 
and focusing on the overdamped regime (rj^ ^ ^(t)) the relaxation time is given by 
T — ?]/S. It follows that the average absolute value of the winding number scales as 



1/8 



W)-l/^f?l . (13) 



In the following section we shall compare this prediction with the numerical solution 
of the set of coupled Langevin equations describing the microscopic dynamics of the 
ion chain. 



4. Results and discussion 

To study the dynamics of defect formation in the structural phase transition we solve 
the N coupled Langevin equations associated with equation ([1]) in the presence of laser 
cooling. This is done in order to obtain reliable predictions for finite systems. We then 
compare the results with the predictions of the continuum field description in equation 
(|4|. whose validity is restricted to the thermodynamic limit. The molecular dynamics 
simulations were carried out using the Protomol platform |32^. Protomol was 
already used before to study the dynamics of cold ion plasmas ^3\. The impulse 
integrator for Langevin dynamics [33] was used. The Coulomb forces were evaluated 
directly between all ions without any approximations. The numerical experiments were 
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Figure 3: Evolution of the order parameter A{x) = \A{x)\e'^^^^^ in a single realization. 
The ion chain is prepared in a linear structural phase of circumference C for which 
the order parameter A vanishes up to thermal fluctuations, i.e. the average off-axis 
displacement of the ions is zero. As the transverse confinement is decreased during 
the quench, the system crosses the critical point and spontaneous symmetry breaking 
occurs, leading initially to the formation of different zigzag domains characterized by 
non-zero \A\ and the choice of an approximately constant value of 9 associated with the 
plane on which the domain lies. The realization corresponds to the exact simulation 
of iV = 400 ions and a quench rate ln(l/a;oTQ) = —5.14, the result is smeared in the 
surface plots to mimic a continuous field, A{x). At longer times (see figure Ufb)) the 
system relaxes to a helicoidal ion chain with non-zero winding number. 




40 aO 12D leO 230 _Q/2 
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c/2 
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Figure 4: (a) Decay of the winding number in a helical chain of 400 ions and 
circumference C. After the transient dynamics associated with the growth and 
formation of the helical chain, the winding number W stabilizes, and at zero 
temperature it remains stable on topological grounds. At finite temperature, thermal 
fluctuations can force the order parameter A{x) to vanish at a given position x during 
the stochastic dynamics. As a result, its phase 9{x) becomes undefined, paving the 
way to the unwinding of the chain. This single realization shows the decay of a = 3 
configuration to a helix with W = 2 which happens shortly after the phase transition, 
(b) Long time evolution of the order parameter 0{x,t) that shows the unwinding of 
the helical chain. 
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Figure 5: (a) Scaling of the winding number W (averaged over ~2000 realizations) 
with the quench rate. The Kibble-Zurek mechanism predicts a power law scaling of W 
with the quench rate, which is quantitatively reproduced by the molecular dynamics 
simulations for moderate quench rates. The linear regression fit gives a power-law 
exponent of 0.129 (regression coefhcient is 0.92) very close to the KZM prediction 
(1/8), while the absolute value {\W\) is about five times smaller than that dictated by 
KZM. At high rates the scaling breaks down and the winding number saturates as a 
result of finite size effects. The length of the vertical lines is related to the standard 
deviation, <jw{tq), of the distribution of W{tq) and is given by hi{aw(TQ))/5, where 
division by five is for the ease of visualization, (b) Normalized histograms for the 
distribution of W corresponding to a quench rate in the parameter region well- 
described by KZM scaling, ln(l/cjo''"Q) = —7.51 (left), and a faster quench rate close 
to saturation region where the KZM scaling breaks down, \n{l / ujqTq) = —3.53 (right). 



done as follows. N = 400 ions are initially set to be equally spaced on the x-axis with 
zero velocity. The simulation is then run for 4x 10'* timesteps, which was sufficient for 
the thermalization of the chain. Then, the transverse frequency is decreased linearly 
between an initial (v^^^) and final (v^^^) value across the critical point j^j'^^ in a time tq. 

This leads to the effective S(t) = SoI/tq with — '2'^t^\i^t^^ "'^t"'')- The simulation is 
continued at the final transverse trapping frequency for another 8x 10^ timesteps which 
ensures relaxation of the ion chain to a configuration in which the winding number 
stays constant. The characteristic frequency of the system is obtained by comparing 
the Coulomb energy and the trapping energy and is given by ujq — \J /AiTeoma^, 
where m is the mass of the ion and a is the spacing between the ions. In terms of ujq 
the parameters used in the simulation were vj.^^ = 2.54a;o, = 1.68a;o, 77 = 4.38a;o 
and ksT — 3.5mujQa^ . 

The evolution of the order parameter A in a typical run is shown in figure [H We 
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observe the formation of correlations near the critical point which leads to regions of 
zigzag chain at constant orientation 9. We can also see points were the displacement 
of ions, \A\, stays zero for considerable time even after the system is quenched through 
the critical point. At these 'phase slip centres' the phase 9 is undefined and so it can 
'slip' by 27r changing the winding number. We have carried out simulations of the long 
time dynamics of the helical chain in order to check whether the winding number stays 
constant. It is observed that the winding number can change shortly after crossing the 
critical point and that it eventually stabilizes. Just after crossing the critical point the 
magnitude of the order parameter A{x) is small enough for the thermal fluctuations 
to bring it to A{x) « at a certain position x. This facilitates phase slips that unwind 
the chain and minimize its Coulomb energy. Figure |4] illustrates the decay of a 14^ = 3 
to a W = 2 configuration. The non- universal nature of this process is beyond the scope 
of the KZM, and deviations between the observed averaging winding number and the 
KZM prediction are expected whenever the decay induced by thermal fluctuations 
plays a relevant role. 

Simulations were carried out for 27 different quench rates. The number of runs 
for each quench rate was ~2000. The results are shown in figured The dispersion of 
the resulting winding number W exhibits a dependence on the quench rate which is 
well captured by the power-law predicted by the KZM. A quantitative agreement is 
found for slow to moderate quench rates i.e. ln(l/a;oTQ) < —4. While KZM captures 
the right dependence on the parameters of the system and the quench, as can be 
seen by comparison of the KZM prediction in equation (fT3|) with the results of the 
simulations, we find that KZM overestimates the winding number by a factor of 5.2. 
This agrees with the fact that the KZM prediction generally surpasses the observed 
density of defects by a factor of up to an order of magnitude, as has been observed in 
a variety of systems [35], including ion chains [5l|6]. Furthermore, the scaling breaks 
down for faster quenches due to a saturation of {\W\) resulting from finite size of the 
chain which leads to a broadening of the distribution. 

5. Conclusions 

Ion chains with periodic boundary conditions undergoing the linear to zigzag structural 
transition in the presence of laser cooling constitute an experimentally attractive 
platform for studying the dynamics of spontaneous symmetry breaking. The critical 
point is well characterized and the transition can be driven by weakening the transverse 
confinement at a given quench rate. We have shown that at an early stage of the 
quench, the resulting chain exhibits multiple zigzag domains oriented along different 
planes in the transverse direction. In a subsequent stage, the topology of the trap 
forces the system to relax to a helical ion chain with a given winding number which 
depends on the quench rate. Furthermore, this dependence agrees quantitatively 
with the power-law predicted by the Kibble-Zurek mechanism with mean-field critical 
exponents for moderate rates. At high quench rates the winding number saturates 
leading to a breakdown of the scaling. This work shows once more that ion crystals are 
very well suited classical simulators of complex and critical dynamics, from condensed 
matter to cosmological models. 
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